DEPARTMENT OF AEROSPACE ENGINEERING 
COLLEGE OF ENGINEERING AND TECHNOLOGY 
OLD DOMINION UNIVERSITY 
NORFOLK, VA 23529 


SUBSONIC AND SUPERSONIC JET NOISE CALCULATIONS 
USING PSE AND DNS 


By 

Dr. P. Balakumar, Principal Investigator 
Farouk Owis 

Department of Aerospace Engineering 


FINAL TECHNICAL REPORT 

For the period ending September 30, 1998 


Prepared for 

NASA Langley Research Center 
Attn.: Ms. Barbara S. Thomson 
Grants Officer 
Mail Stop 126 
Hampton, VA 23681-2199 


Dr. Michele Macaraeg, Technical Officer 


Under 

Research Grant NAG- 1-2054 
ODURF File No. 182431 


March 1999 



Jet stability and Noise Prediction 
Using Parabolized Stability Equations 


i 



Prediction of Supersonic Jet Noise 


P. Balakumar 

Aerospace Engineering Department, 

Old Dominion University, Norfolk, Virginia 23529-0247 


Abstract 

Noise radiated from a supersonic jet is computed using the Parabolized Stability Equations 
(PSE) method. The evolution of the instability waves inside the jet is computed using t e 
method and the noise radiated to the far field from these waves is calculated by solving the 
wave equation using the Fourier transform method. We performed the computations for a cold 
supersonic jet of Mach number 2.1 which is excited by disturbances with Strouhal numbers St-2 
and 4 and the azimuthal wavenumber m=l. Good agreement in the sound pressure level are 
observed between the computed and the measured (Troutt and McLaughlin 1980) results. 


1. Introduction 

In this work we computed the noise radiated from supersonic turbulent jets using the PSE 
(Parabolized Stability Equations) method. Jet noise can be divided into three categories: ( ) 
shock-induced screech tone noise; (2) shock-induced broad band noise; and (3) turbulent mixing 
noise The screech tone appears in an imperfectly expanded jet as discrete band in the front part 
of the jet. The screech tone phenomena is very complex and it is believed that the toroidal and 
helical vortices which shed at the lip of the nozzle interact with the shocks when they propagate 
downstream and makes the shock to oscillate and this radiates sound in the upstream direction 
at discrete frequency. In an imperfectly expanded jet, the shock-cells formed by the oblique 
shocks or the expansion fans generated at nozzle lip interact with the large scale turbulence 
and generate broad band noise. The dominant part of the broad band shock-associated noise 
is comprised of a spectral peak with a relatively narrow half-width. Turbulent mixing noise is 
the noise component which is contributed from the large-scale and small scale turbulence in the 
jet For a perfectly expanded supersonic jet, the noise is completely generated by the turbulence 
in the jet and the predominant part of the noise is radiated in the downstream direction in 
the range between 25^5°. It is observed that at low Reynolds numbers the noise is radiated 
at a discrete frequency which is closer to the most unstable instability wave for that jet. A 
moderate and high Reynolds numbers there is discernible peak but they become broad band 
These similarities between the noise generated from the high and low eyno s num er J 
imply that the noise generation mechanisms in supersonic jets are same at different Reyno 
numbers. Several experiments were performed to identify these mechanisms ( McLaughlin e a • 
1975 1977, Morrision and McLaughlin 1979, 1980, Troutt and McLaughlin 1982 1, Seiner et al. 
1993) and it is concluded that the dominant part of the turbulent mixing no.se of high Reynolds 
number supersonic jets is generated by the large-scale coherent structures. It is also concluded 
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that these large-scale coherent structures are the instability waves of the jet (Tam 1971, 1972, 
Morris and Tam 1977, Tam et. al. 1991). 

The jet column can be divided into three regions, namely, laminar, transitional and turbulent. 
The extent of the laminar and transitional region depends on the Reynolds number and the Mach 
number. The laminar region extends a few diameters from the exit and the disturbances grow 
exponentially in this region. Further downstream, when the amplitude of these disturbances 
becomes large, they interact nonlinearly and then saturate. These instability waves then disengage 
their coherence and disintegrate, and the flow becomes turbulent. In a senes of detailed 
experiments (McLaughlin et al. 1975, Morrison and McLaughlin 1979, 1980, Troutt and 
McLaughlin 1982), it was confirmed that most of the noise is generated in this region where 
the disturbances saturate and decay. It is also identified that this region is always close to the 
end of the potential core region. 

Several methods are explored in practice to compute the turbulent mixing noise generated 
from supersonic jets. One is the multiple scale approach ( Tam 1991). In this method it is 
assumed that the flow field induced by the large scale structures in a turbulent jet can be modelled 
as the flow field generated by the evolution of instability waves in a given turbulent flow. The 
mean flow is obtained from the experimental measurements or by solving the Reynolds averaged 
Navier-Stokes equations. The local growth rate, the wavenumber and the eigenfunctions of a 
disturbance with a constant frequency and a fixed azimuthal wavenumber are obtained by solving 
the compressible Rayleigh’s stability equations and the nonparallel corrections are detenmne 
using the multiple scale and the adjoint method. After computing the flow field in the near field of 
the jet the acoustic field in the outer part of the jet is calculated by solving the wave equation. 
In the second method, the mean flow field and the flow field associated with the large sea e 
structures are computed by solving the Large-Eddy-Simulation equations with the subgrid scale 
modelling (Mankbadi et. al. 1994). The noise in the far field is again obtained by solving the 
linearized Euler equations. In the third method, the complete Navier-Stokes equations are solved 
(Mitchell et. al. 1996) to obtain the inner and the outer acoustic fields. The last two methods are 
more accurate and do not make any assumptions other than the subgrid scale modelling in t e 
LES approach. However, they are computationally very expensive compared to the first metho . 

In this work we followed the first approach and compute the flow field in the inner part of the 
jet using the PSE method instead of using the multiple scale method and compute the acoustic 
field in the outer region by solving the wave equation with the pressure obtained from the P 
as the inner boundary data. We obtained the mean flow by curve fitting the experimental data 

by cubic splines. 


2. Formulation 

We are concerned with the evolution of a small disturbance of a single frequency and a fixed 
azimuthal wavenumber inside an axisymmetric supersonic jet and the noise radiated from this 
disturbance field to the outer part of the jet. Figure 1 shows the schematic and the coordinate 
system that we used in the analysis. We perform the computation m two steps. In the first step, 
we compute the flow field inside the jet using the PSE method and in the second step we solve 
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the wave equation to compute the acoustic field in the far field of the jet. Since we are seeking 
linear solutions, the disturbance quantities can be written in the form 


im6 — iut 

q(x,r,6,t) = q(x,r)e 


( 1 ) 


where x, r and 0 are the axial, radial and the azimuthal coordinates; m is the mode number in 
the azimuthal direction, u is the frequency and t is the time; q is a vector, 

q = {u,t>, w,p,T} T , ( 2 ) 

and u, v and w are velocity components in the x, r, and 6 directions, p is the pressure and T is 
the temperature. We non-dimensionalize the variables as follows: 

velocity * Uj : jet exit velocity 

Temperature ” T 0 : ambient temperature 

Density " po : ambient density. ^ 

Mach number based on the jet exit conditions Mj=— ^^=- 

Mach number based on the jet exit velocity and the free stream temperature M o= ^ ^ ■ 


2.1 The Parabolized Stability Equations 

In the parabolized stability equations (PSE) approach, one attempts to construct an approxi- 
mate solution of the full Navier-Stokes equations. The idea was first introduced by Herbert(1991) 
and applied to linear and non-linear Blasius boundary-layer flow by Bertolotti (1992). Now it 
has been developed and has been applied to two and three-dimensional incompressible an 
compressible boundary-layer flows and supersonic jets (Chang et.al. 1994 , Malik et.al 1994 
Balakumar 1994, Malik et.al 1997). We give a brief description of PSE for a general flow and 
apply to jets as a special case. 

In the linear parallel stability analysis, the disturbance quantities are written in normal mode 
form. If q(x,y,z,t) is a flow variable in normal mode analysis, we write 


i f adx + iflz — iuit 

q(x,y,z,t) = q{y)e 


( 3 ) 


Here u, is the frequency of the disturbance, f3 is the azimuthal or spanwise wave number, and 
o(x) is the wavenumber in the axial direction. The eigenfunction is q( y) which is function 
of y (normal coordinate) only. We substitute this expression into the linearized Navier-Stokes 
equations and, assuming the flow is parallel in the streamwise direction , we obtain an ordinary 
differential equation for q( y). This equation, along with the homogeneous boundary condition 
at the wall (for boundary-layer flow) and in the free-stream, forms the eigenvalue problem for 
the wave number a and for the eigenfunction q{ y). 
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In practice, the mean flow is not constant, but varies in the streamwise direction. Since this 
variation is relatively slow, the wavenumber a and the eigenfunction q{ y) vary slowly in the 
streamwise direction. In the PSE approach, we use this fact to construct an approximate solution 
to the Navier-Stokes equations. In the PSE formulations, we write the disturbance quantities as 


/ ,, w \ t / adx+idz—iwt 

q(x,y,z,t) = q(x,y)e J 


where q(x, y) 
The first 


represents the amplitude part and the exponent as the wave 
and the second derivatives of q(x,y,z,t) are 


part. 


dq_ 

dx 


S^ia{x)q{x,y) + 


8q 

dx 



adx+ifiz 


iu>t 



dx 2 


i / otdx+ipz—iut 
C J 


( 4 ) 


( 5 ) 


(6) 


Since the amplitude part q varies slowly in the x direction in the PSE approximation, we neglect 
and write the second derivative as 


d 2 q 

Yx 2 


o . da 

-a“{x)q{x, y) + V) + 2 ta 

J adx-\-if)z — iuit 


( 7 ) 


Therefore, if we substitute this approximation into the Navier-Stokes equations, the second 
derivative in the x direction drops out and the system of equations is parabolized. The solution 
may be found by marching in the x-direction which is the major advantage of PSE. 

The first step of the procedure is to start, as in any other parabolic problem, with a known 
solution at x=xo, <?(x=xo,y), a(x=xo), assume a(x=x 0 +Ax)=a(x=xo), march the PSE equation to 
the next station x=xo+Ax, and solve for <j(x=xo+Ax,y). The second step is to compute the new 
a(x=xo+Ax) from the computed <j(x=xo+Ax,y). The problem is that since q is a function of y, 
and q can be any physical quantity (e.g., velocity, temperature, mass flow rate, etc.), it is not 
clear how to compute a. This problem arises in the experiments too. If we want to measure 
growth rate using hot wires, we face the same problem as to which quantity and at what location 
we want to measure. Therefore in a nonparallel flow, there is no unique wavenumber as in the 
parallel linear theory. We call the wave number we compute the “wave number based on some 
quantity (e.g., velocity, temperature, mass-flow rate)”. Usually the wave number is computed 
at the location in the shear layer where the flow quantities become maximum. For example, if 
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we choose to compute the wave humber a a. the location where streamwise velocity u peaks, 
a is obtained from the relation 


a(x — Jo + Ax) 


: 


lUmax 


du 

dxj y = y 


( 8 ) 


1 rz T In the third step after the newa(x=x 0 +Ax) 

where y is the location where u pea s . . are repeat ed at this location with the 

“ SS de'cribe the goveming .natrons and 

- — -r — - - -* 

Let X|, x 2 , X, be a set of generalized ortnogon The mean 

coefficients are h„ h 2 , h, and the velocmes ,n the threchon quantitie s are 

how quantities are represented *<>-<£' vanablcs are 


flow quanut.es are . //.Therefore, the total variables are 

represented by Q — \U\^ u 2'> u i‘> a ’ 5 


Q = Q° + Q' 


(9) 




, - . i f adx\ (10) 

Q> = Q(xi,x 3 )e J 

substituting this expression into the goveming equations for <? and neglecting the nonlinear 
terms we obtain the following linear PSE equations for Q, 


A d S + B^ + CQ = 
dx\ ox 3 

D dx\ + dxidx 3 


(ID 


Here, A, B, C, D and IE are I ^".liS "erW^ by 

two ^oim^pwtnd scheme 3 amf'second by applying fl^variableTand 6 the^metrics 

the norma, direction When we apply Ihrs to an axr.ymmetrrc ^ ^ ^ az|muthal 

becomes xi=x, X 2 -$> x 3- r > “i * 2 ’ 3 

and in the radial directions are u, w and w appropriate boundary conditions and 

the governing equa^n^h^^to^be derived^separately. We use L’ Hospital’s rule at r=0 and obtain 
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T'T and ,he b0Undary COndi,ionS The bound ^ conditions at 
different form based on the azimuthal mode number m. They take the following form: 


m = 0 




du 

a7 = °' 

II 

o 

S 

II 

o 

86 

(12) 

m = ±1 




u = 0, 

• rn 

W + 1 — '-V = 0, 

m 

9 = 0. 

(13) 


rn > 2 

u = 0, v = 0, w = 0, 9 = 0. (14) 

In the far field r — ► oo, we use the condition 


u = 0, v = 0, w = 0, 9 = 0. 

(15) 

of the ™ mpUtati ° nS ’ WC ' mpOSe these boundar y conditions at r=50dj, where dj is the diameter 


2.2 Wave Equation 


the 

the 


To compute the acoustic field in the far held, we solve the linearized Euler equations using 

Fourier Transform technique. In cylindrical coordinates the linearized Euler equations for 
amplitude part q(x, r ) become 


. - 1 dp 

■ IUJV = - 

p dr 
. _ im 

ILUW = — — — O, 

pr 

. . 1 dp 

turn = — , 

p dx 

~ , 1 ( v , dv im du 

lUP+ w{r + lTr+T W+ Tr 


(16) 


= 0. 


Eliminating the velocity components , we obtain the wave equation for the pressure 


d 2 p d 2 p 1 dp 

d? + M + + l pM ‘ u 


p = 0. 


(17) 
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The boundary conditions are 


p — » 0, as r — > oo, 
p(x,r 0 ) = po(x), at r = r 0 . 

Here ro is the height at which the lower boundary for the wave equation is located and this line 
is in the region where the mean velocity becomes almost zero ( Figure 1). Usually, we take ro~ 
2-3 diameters of the jet. po(x) is the pressure distribution along the line r=r 0 which is known 
from the PSE calculations. The wave equation (eq. 13) is solved using the Fourier transform 
technique and the solution is given by 



Po(k x ) 


X m{ T ^ ) 
k-mij' 0 ^ ) 




dk r 


A 2 = k 2 x - pM 2 u> 2 , (Real(X) > 0). 


(19) 


Here K m is the modified Bessel function of second kind of order m and the integration is along 
the contour C, which does not cross any branch cuts in the complex k x - plane and is shown in 
Figure 2. We evaluated this integral numerically. For comparison and to determine analytically 
the intense noise radiation direction and the dominant wavenumber region, it is advantageous 
to evaluate this integral asymptotically for large distances. As it was done by Tam and Burton 
(1984), we rewrite this expression in spherical coordinates (R, X, <f>). In the new coordinate 
system, the pressure expression becomes 



Pojkx) 

X m ( T 0 X ) 


K m (RsinX X)e' kxCOsX R dk x . 


( 20 ) 


For large R, the Bessel function and the integral become 


,< mi Rs,nXX) - 
P{x,r) = 


2 -RsinXX 


/ K [ Po(kx) 1 i{k x cosX+tsmX\)R ,, 

V 2 RsinX J K m (r 0 A) 


( 21 ) 


This integral can be evaluated using the stationary phase method and the integration yields 

P (R,xA,i) = U- tii ; M “ c ° aX) 1 

R \ K m {r§MujsinX) f (22) 

MuR+m(f)— wt— 
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and the stationary point is given by 


k x — MuicosX. (23) 

From this expression for a fixed R we can determine the maximum acoustic pressure and the 
direction of the maximum pressure and the wavenumber by evaluating the expression in the 
bracket in eq. (22).. 


3. Results 

Computations are performed for a cold supersonic jet of Mach number Mj=2.1. The mean 
flow velocity profiles are obtained by curve fitting the experimental data which is given in 
McLaughlin, Seiner & Liu (1980). We used the same empirical formulas as used by Tam & 
Burton (1982). The jet is divided into three regions as core, transitional and fully developed 
regions. In each region, different functions are used to represent the measured velocity profiles. 
The functions used are: 


Core region : 0 < x < 


U = 



(r < h), | 
{h < r). J 


(24) 


Transition region : 

Xf < X < X f 

_ ( u c {x) 


1 u c (x)exp 

[-«» 2 )(^) 2 ] 


(r < h),} 
(h < r). J 


(25) 


Fully developed region : x > x f 
u = {u c (;r)exp - ( /n2 )(i(iy) 1 }• 


(26) 


Here h(x) is the thickness of the core region, b(x) is the height from the end of the potential core 
to the half velocity point, U=.5. x t and Xf are the locations of the end of the core and the start of 
the fully developed regions from the nozzle exit. From the experimental results of McLaughlin, 
Seiner & Liu (1980), it is found that the core region extends up to five diameters and the fully 
developed region starts at eight diameters. b(x) is determined by approximating the measured 
thickness of the shear layer by a cubic spline curve. The variation of h(x) in the core region and 
the variation of u c (x) in the fully developed region are obtained from the conservation of axial 
momentum condition. In the transition region, the variation of h(x) and u c (x) are obtained by 
matching their absolute values and their derivatives at both ends. Figure 3 shows the measured 
and the curve fitted distribution of b(x) and figure 4 shows the corresponding variation of the 
centerline velocity u c and the thickness of the core h(x) along the axis of the jet. The temperature 
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distribution is obtained by assuming that the total temperature is constant across the jet and the 
radial velocity is obtained by integrating the continuity equation in the radial direction. 

After the mean flow profiles are obtained, the evolution of the disturbances inside the jet is 
computed using the PSE method. In all the previous computations, Euler equations are solved 
to obtain the flow field inside the jet. The main reason is that we are simulating a turbulent jet 
and it is difficult to define a Reynolds number in this case and further that the flow dynamics 
are of the inviscid type in free shear flows and the Reynolds number has little effect in the 
results. The PSE equations are derived from the complete viscous Navier-Stokes equations and 
we performed the computations as a viscous problem but fixing the Reynolds number arbitrarily 
at larger values ~ 10 6 . We present the results for two Strouhal numbers St=.2 and .4 and for the 
helical mode m=l which is the most amplified disturbance in an axisymmetric jet. 

Figure 5 shows the evolution of the mass velocity fluctuations with the axial distance obtained 
from the nonparallel PSE calculations for the Strouhal numbers St=.2, .4 and m=l. The mass 
velocity fluctuations are measured at the radial location where they are maximum. Figures 6 
and 7 show the real part and the amplitude of the pressure distribution in the axial direction at 
the radial location j- = 3 for two Strouhal numbers St=.2 and .4. We observe that the pressure 

Y * 

peaks at about j- = 8 and decreases further downstream and we also notice that the wavelengths 
are smaller at St=.4 than that at St=.2. 

Figures 8 and 9 show the spectral distribution of the pressure for the Strouhal numbers 
St=.2 and .4 obtained by taking the Fourier transform of the pressure distribution shown in the 
figures 6 and 7. The branch points in the complex k x plane are located at k x = ±u>M. For 
St=.2 and .4 these points are located at k x =± 1.922 and 3.844 respectively. The waves with the 
wavenumbers larger than these values travel subsonically relative to the free stream and hence 
will not radiate noise to the far field. Equation (22) shows that the far field noise is determined 
by the quantity and the noise radiation direction and the wavenumber is related by the 

eq. (23). In figures 10 and 1 1 we plotted these two quantities for the Strouhal numbers St=.2 

O 

and .4 respectively. The figures show that the intense noise is radiated at an angle of 45 from 

0 

the axis for the frequency St=.2 and the angle is 37 for St=.4. Figures 12 and 13 show the 
sound pressure levels (SPL in dB) radiated from the instability waves with the frequency St=.2, 
.4 and m=l to the far field obtained by integrating the eq. (20) numerically. We see that the 
noise is radiated in a wedge shaped region and the intense noise is radiated in a fixed direction 

O O 

which is inclined at an angle of 44 for St=.2 and is inclined at 37 for St=.4 which agree with 
our earlier prediction from the asymptotic theory. Figures 14 and 15 show the experimental 
measurements of sound pressure levels by Trout and McLaughlin (1982) for a jet excited at 
St=.2 and .4. In figures 12 and 13 we matched the computed sound pressure levels with the 
experimental measurements shown in figures 14 and 15 at one point = 30 and -fi- = 20. 
It is seen that the agreement between the experiment and the computation is very good. In 

O 

the experiment the noise is radiated at an angle of 44 and 36 which are exactly same as that 
predicted from the computations. To understand the character of the flow field inside the jet, in 
figure 16 we plotted the variation of the quantity 



along the axis. Here c is the phase speed of the instability wave and a<x> is the speed of sound 
in the freestream. If this quantity is positive the wave is travelling subsonically relative to the 
freestream and if it is negative it is travelling supersonically. We observe from the figure that 
the waves travel supersonically for the first 5 and 7 diameters for St=.2 and .4 respectively. 
After that they travel subsonically and further downstream they again travel supersonically. This 
figure depicts that the most of the noise is radiated from the first few diameters of the jet 


4. Conclusions 

In this work, the evolution of a small disturbance of a single frequency and a fixed azimuthal 
wavenumber inside an axisymmetric supersonic turbulent jet is computed using the PSE method 
and the noise radiated from this disturbance field to the outer part of the jet is evaluated by 
solving the wave equation using the Fourier transform method. The PSE computations are very 
efficient and take about 5 minutes on a SUN workstation. Computations are performed for a jet 
of Mach number 2.1 at two Strouhal numbers St=.2 and .4. Good agreements are found between 
the computed and the measured sound pressure levels. It is also observed that the waves travel 
supersonically for the first few diameters from the jet exit and after that they travel subsonically. 
This agrees with the observation that most of the noise in a supersonic jet is radiated from the 
first few diameters of the jet. 
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Figure 3. Axial distribution of mean velocity profile 
parameters b(x). 



Figure 4. Axial distribution of mean velocity profile 
parameters b(x), h(x) and U c - 
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Figure 5. Axial distribution of mass velocity fluctuations in 
the shear layer for St=.2 and .4 and m=l. 



Figure 6. Pressure 


t=.2 and m=l. 
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Figure 9. Spectral distribution of pressure p 0 (k*) for St-4 and m-1 
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Figure 11. Distribution of po(k,)/K m0 (k,) and the sound radiation 
angle given by eq.(23) for St= 4 and m=l. 



Figure 12. Sound pressure level contours computed from the wave 
equation for St= 2 and m— 1. 
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Figure 15. Sound pressure level contours measured by Troutt and 
McLaughlin (1982) for St=.4 and m=l. 



Figure 16. Axial distribution of the quantity (k*-u;M) obtained from 
PSE for St=.2 and .4. 




Jet stability and Noise Prediction 
Using Direct Numerical Simulation 


i 



Contents 


1. Introduction 

2. Governing equations 

3. Numerical schemes 

4. Boundary conditions 

5. Noise calculations 

a) Mathematical formulations 

b) Numerical discretization 

6. Results 
References 


2 



1. Introduction: 

The development of high speed civil transport plane requires reducing the jet exhaust 
noise. The computations of the jet noise radiated from the mixing layer may be divided 
into two parts. The first part is the evaluation of the near-field source and the second part 
is the computations of noise in the far field. In this work, we concentrate on obtaining the 
sound source using the direct numerical simulation of the full unsteady Navier-stokes 
equations. 

Experiments by Laurence (1996) have shown that the sound power emitted from the jet is 
greatest within 4 to 5 diameters downstream and then decays through a transition region. 
This region is characterized by large vortical structures and is not fully turbulent (Soh 
1994) which gives the motivation that we solve the unsteady flow equations using the 
direct numerical simulation (DNS) to provide the sound source for an acoustic 
computation of the far-field noise. Due to the limitations of the computational facilities 
available at the present time, the linearized Euler equations or the linearized wave 
equation is used to calculate the noise in the far field. The linearized Euler equations 
approach neglects both viscosity and nonlinear effects. The viscosity effects can be 
neglected since the free shear layer in the far field is essentially inviscid (Mankbadi 
1992). The nonlinear flow effects and source generation are confined to the near field 
(Shih 1995) and sound propagation in the far-field can be modeled by linear instability 
waves. Tam and Morris (1980) calculated the noise in the far-field using the linear 
instability waves. 

In order to evaluate the sound source using DNS, the simulation must be performed using 
numerical techniques with minimum distortion and diffusive characteristics. The 
numerical errors get worth for high Reynolds number flow simulations. Typically, free 
shear layer flows of interest have very high Reynolds numbers. Therefore, higher order 
accurate numerical schemes with minimum dissipation and dispersion errors are needed. 
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Part of this study is dedicated to the investigation of the accuracy of different numerical 
schemes. 

The treatment of the boundaries is very important in getting an accurate solution of the 
Navier-stokes equations. Various computational techniques have been developed in the 
past to minimize the reflections of the out-going waves. Some of these techniques are 
based on the characteristics of the equations such as Thompson (1987) and Giles (1990). 
Other methods are based on the far-field asymptotic solution (Bayliss 1980, Enquist 
1979, Hagstrom 1988, Tam 1993 and Tam 1995). In addition, a buffer zone technique has 
been developed in which the mean flow is altered gradually to be supersonic in a buffer 
region adjacent to the artificial boundaries (Givoli 1991, Street 1989 and Ta’asan 1995). 
A different approach has been proposed by Hu (1996) to damp the disturbance 
exponentially in a short layer called perfectly matching layer. We devoted part of this 
work to investigate different types of boundary conditions which will give minimum 
wave reflections near the boundaries. 

In order to evaluate the jet stability and its radiated sound, the jet mean flow has to be 
computed. We used the boundary layer equations to solve for the jet mean flow. The 
equations are solved using two-point compact finite difference scheme. In addition, the 
linear stability theory is used to calculate the most unstable frequency for the jet and its 
eigenfunctions which will be used to excite the jet inflow boundary in DNS code, 
linearized Navier-stokes equations are used for the linear stability analysis. 

It is widely believed that the large vortical structure and the vortex pairings are 
responsible for jet noise radiation (Colonius and et al. 1997 and Mitchel 1995), we 
investigated the sound generated by vortex roll-up and pairings by forcing the inflow 
with the most unstable frequency f and its subharmonics f/2 and f/4. 
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Spectrum 



Figure 7. Pressure distribution at r/dj=3 for St-.4 and m-l. 



Figure 8. Spectral distribution of pressure po(k*) 
for St= 2 and m=l . 




1. Governing Equations: 

Navier-Stokes equations will be used to evaluate the near field source of jet noise and to 
compute the linear and nonlinear instability waves of the mixing layer. The equations are 
written in conservative form, cylindrical coordinates and for 3-d axisymmetric jet. The jet 
radius, exit velocity, temperature and density (q , Uj , Tj ,p ] ) are used to normalize the 
equations. 
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Here Q is the unknown vector, F, G, and H are the fluxes in x, r, 0 directions, 
respectively; S is the source term. The shear stresses are calculated as follows: 
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2. Numerical Schemes: 

MacCormack type schemes with operator splitting and high order accuracy developed by 
Hixon (1997), up to sixth order accurate in space and fourth order accurate in time, will 
be used in this simulation. The operators are applied in the following symmetric way: 

Q n+2 = L lx L 2 r L 2 g L xe L yr L {x Q n 

Where L x , L r and Lq, are one-dimensional operators in x, r and 0 directions that are 
applied to these one-dimensional equations 
Qt = - F x 
Q t = - G r + S 
Qt = - He 

The operators will be alternated with symmetrical variants such that the scheme will 
maintain the accuracy. Let Li be the operator with forward finite difference in the first 
step, then L 2 is the operator with backward finite difference in the first step. 

Using Rung-Kutta, the equations are integrated in time as follows: 

o,=~[Fm 

ax 

h,=-At^ [F{Q n )\ 

ax 

j b 

h 2 = -At — [F(Q n + a 2 h^)] 
ax 

h,=-At^- f [F(Q n +a,h 2 )] 
ax 

h 4 = -At [F((2 n +«A)] 
ax 

= -At -y-' [F(Q n +a 5 h 4 )] 
ax 

j b 

K ~ ~At — [F(Q n +a 6 h 5 )] 
ax 

Q n+ ' =Q n + 0^+ 0 2 h 2 + 0^ + 0 A h 4 + 0 5 h 5 + 0 6 h 6 
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Accuracy 

a 2 

OC 3 

a 4 

a 5 

0t6 

Pi 

P 2 

p3 

p4 

Ps 

P6 

Second 

order 

l 

0 

0 

0 

0 

1/2 

0 

0 

0 

0 

0 

Fourth 

order 

1/2 

1/2 

1 

0 

0 

1/6 

1/3 

1/3 

1/6 

0 

0 

Multi-step 

(4-6) 

First step 

1/2 

1/2 

— 

0 

0 

1/6 

1/3 

1/3 

1/6 

0 

0 

Second 

step 

.3533 

.9996 

.1522 

.5342 

.6039 

.0468 

.1373 

.1710 

.1976 

.2823 

.1651 


Space discretization: 

MacCormack schemes with higher order accuracy up to sixth order and minimum 
dispersion error are used. The forward and backward discretization for the fluxes are 

written in this form. 


j = a _, F w + a 0 F, + a, F M + a 2 F i+2 + a 3 F M 
= a -\F M + a 0 Fi + a,F M + a 2 F „ 2 + a,F t _, 

\ dx j 

The accuracy of MacCormack schemes is increased to sixth order accurate by adding a 
point ( a 3 ) to the discretization in the backward and forward directions. Also, the 
dispersion error is optimized for these schemes by adding another point (a _,) such as 4/4, 


6/4 and DRP/opt. Schemes as shown in the following table. 


Scheme 

a .1 

ao 

ai 

a 2 

a 3 

" 2/2 

0 

-1 

1 

0 

0 



Ax 

Ax 



4/2 

0 

-7 

8 

-1 
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6Ax 

6Ax 

6Ax 


6/2 

0 

-37 

45 

-9 

1 



30 Ax 

30Ax 

30Ax 

30Ax 

4/4 1 

-2 
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6Ax 

6Ax 

6Ax 


6/4 

-9 

-19 

36 

-9 

1 


30Ax 

30Ax 

30Ax 

30 Ax 

30Ax 

DRP/opt 

- .30874 
Ax 

-.63254 

Ax 

1.2330 

Ax 

1 -.3334 

Ax 

.04168 

Ax 
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3. Boundary conditions: 

4.1 Inflow B.CS: 

For supersonic flow, all the flow variables are specified at the inflow boundary using 

Q = Q m +£ Qdist 

Jt=0 

where 

Q m is the mean flow variables, 

Q k ( r ) are the eigenfunctions of the flow variables which are obtained from the 
linear stability code for the most unstable frequency to and its subharmonics to/2 and co/4 
and n is the azimuthal wave number, e is the amplitude of the disturbance at the inflow 
boundary. 

For subsonic flow, one flow variable is obtained from the interior points using the 
continuity equation and the rest of the flow variables are specified at the boundary as 
specified above. 

dp dp u + dp v + 1 9/9 w _ p v 
dt dx dr r 96 r 


The continuity equation may be written in the following characteristic form: 


3P« +I ^ + i 


where 


9r 


dx dx 

2 dp dp 
A 2 — u (c ) 

dx dx 

, dp 9 u „ 

A 5 = (m + c) (— + pc — ) 
9jc dx 


dr 


96 


For non-reflecting inflow boundary: 

\ 2 =\ 5 =0 and X\ is calculated from the interior points. 
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4.2 Centerline Conditions: 

Different techniques will be applied to overcome the singularity at the centerline for 3-D 
jet: 

• Using the average values of the flow variables on the ring r=Ar. Prior to taking the 
average on the ring r=Ar, some conditions have to imposed on the ring r=0 based on the 
physics of the flow at the centerline. The first condition is obtained by noticing that the 
multivalued nature of the flowfield doesn’t extend to pressure, temperature, density and 
axial component of the velocity. The radial and azimuthal components of velocity have 
multivalued at the centerline ( Griffin et al. 1979). The second condition derives from the 
fact that the true velocity vector through any point can have only one direction in 
physical space, Since in our problem there is a symmetry plane (x, r plane) through which 
there can be no mass flow. If the velocity vector at a given point x on the centerline has 
the value U avg , then the radial and azimuthal components (v , w) must be: 
v(x,r,6) = U avg cos 6 
w{x,r,6) = ~U avg sin # 

To determine U avg , one computes the average of the velocity vector U 
On the ring using 

U(x,r,6) = v(jc,/\6)cos6 -w(x,r,6) sm6 
The average on the ring r=Ar is computed using this equation: 

j k max 

q ,_ w = Y q , j=2, 1=1 ,imax and 1= 1 ,kmax 

’ fcmax n ijx 


where q = (puU pT E) 

• Considering the centerline as an interior point with r=e and using the same governing 
equations. 

• Using Navier-stokes equations in Cartesian coordinates at the centerline only and using 
the following transformations to switch from Cartesian to polar coordinates and vice 
versa. 


L U z J L- 


cos# 

sin# 


-sin# 

cos# 


cos# 

-sin# 


sin# 

cos# 


10 



For 2-D jet, a new set of equations will be derived using L’Hospital rule. These 
equations will be applied with the following symmetric conditions: 

dp _du _ dp _ 
dr dr dr 

v = 0 


4.3 Outflow and Radiation B.CS: 

Different outflow and radiation B.CS will be investigated in this study. 

i. Thompson’s Characteristic B.CS: 

The governing equations (1) could be written at x=L in the following characteristic 
form: 


dQ , 1 drG 1 dH _ _ 

^ +rf+ + = c j+ s 

dt r d r r do 

where d is the amplitude of the characteristic waves and C is the shear stresses in x- 
direction. 


dj = [A^ + (A, + Aj ) / 2 ] / C i 2 
d 2 = wd, +(A 5 -A,)/2C 

d 3 = vd { + p A 3 

d 4 = wd, + p A 4 

d 5 =-^(u 2 + v 2 + w 2 )*d t + p ud 2 + p vd 3 + p wd 4 + (A, + A,)/[2(/-l)] 


where 


- . .dp du 

A, =(u-c ) (- — pc—) 

ox ox 

^ =U Y X 


, dp dp, 


* 2=u(c ih-£ ) 




dw 



i . dp du 

A 5 =(u + c ) (- — pc—) 
ox ox 

For nonreflecting outflow B.CS, all the incoming characteristics are set equal to zero 
and the outgoing characteristics are calculated from the interior points. Similar 
equations are obtained for the characteristics running in r-direction at r ma x- In the 
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plane of intersection, both characteristic equations are applied. In the plane of 
symmetry at 0=0 and 0=rc, symmetric B.CS are applied. 


ii. Outflow Boundary [proposed by Rudy and Strikwerda (1980)]: 

In these boundary conditions, a pressure correction term is added to the incoming 
waves in order to keep the pressure at the outflow close to the mean value. So, the 
incoming wave will be corrected as follows: 

A 5 = K(p-p m ) 


where 


K = 


g(l~M L) c 
L 


M max is the maximum Mach number in the flow field, C is the local speed of sound, L 
is a characteristic dimension of the domain, and a is nonreflecting parameter ranges 

between zero and one. 


iii. Buffer Domain B.CS: 

The buffer - domain technique was proposed by Streett and Macaraeg (1989). The 
technique is based on gradually reducing the ellipticity of the Navier-Stokes equation. 
The sources of the ellipticity in the equations are the streamwise shear stresses and the 
pressure terms. To deal with these sources, the streamwise viscous terms and the 
pressure derivative in the streamwise direction are smoothly reduced to zero through 
multiplication by the following attenuation function: 

Sj =i[l + tanh(4(.-2^-)ll 

where Nb marks the beginning of the buffer domain and N x marks the outflow 
boundary location. 


12 



iv. Perfectly Matching Layer: 

In this technique, a region is added at the boundary as shown in figure (4.1) to damp 
the disturbance. In this region of the domain, exponential damping terms are added to 
the governing equation (2.1) of the form: 


M + *Z.T ^> +g , gl=0 

dr dx 

dQ2 , 1 dr(G-G m ) 


dt r dr 

Q = Ql + Q2 + Q„ 


+ <r x Q2 = (S-S m ) 


where 



b. 


where a mx is the maximum value of a, Lb is the length of PML and x* is the beginning 
of the Layer. Q„,F m , G m and S m denote the mean flow variables. 


a x =0 


a y*0 

a y*0 

a x =0 

<7x^0 

Q 

n 

o 

Gy=0 


Figure (4.1): Computational Domain 
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5. Noise Calculations: 


a. Mathematical formulations: 

The jet noise and the pressure disturbance will be calculated in the far field using the 
linearized wave equation. The linearized Euler equations in cylindrical coordinates may 
be written as: 


du _ 1 dp 

dt p m dx 
9v _ 1 dp 

dt P m dr 
dw _ 1 dp 

dt rp m de 

dp 1 , 1 9w v dv du . n 

' r(-~ +-+— + — > = 0 


dt + M 


r dO r dr dx 


(5.1) 


Eliminating the velocity components from the last equation, we obtain the wave 


equation 

d\p 1 , 1 d\p 1 9v a 2 /? d^_p 

dt 2 p m M 2 o r 2 d0 2 rdr dr 2 dx 2 


(5.2) 


where p m is the free stream mean density 


u . 

and M a - = ,T 0 is the free stream temperature 

The wave equation (5.2) is solved with boundary conditions at r=r 0 and at the far field. 
At r=r 0 the pressure is known as function of time, azimuthal direction and axial position 
from the direct numerical simulation code and the pressure is applied as Dirichlet 
boundary condition. 

p 0 = p(x, r o ,0, t) (5-3) 

At all other boundaries, the following radiation boundary condition derived by Tam 
(1980) is applied: 

JJp + 3p +z = 0 (5.4) 

C a dt dR R 

Where R is in the spherical coordinates and C a is speed of sound in the far field. 
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R = y x 2 +r 2 

Then equation (5.4) may be written in cylindrical coordinates as follows: 

_L3p + £^p + ^ap + £ =0 ( 5 , 5 ) 

C„ dt R dx R dr R 

For two-dimensional axisymmetric jet, the inhomogenous term is divided 
by 2R rather than R (see two- and three-dimensional conditions by Roe 1989). 

b. Numerical discretization: 

The wave equation (5.2) is solved along with the boundary conditions (5.3) and (5.5) 
using six order compact finite difference Pade scheme and fourth order Runge Kutta is 
used to integrate the derivatives in time. 
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6. Results: 

The results are presented here for two-dimensional subsonic and supersonic 
axisymmetric jet with high and low Reynolds numbers. Due to the limitations of the 
computational facilities, three-dimensional jet simulation is under development. The 
simulation is done for high subsonic jets at Mach numbers 0.8 and 0.85 while the 
supersonic jet results is presented for Mach number equal 2.5. 

Figures (1) through (3) represent the variation of the mean flow axial velocity, radial 
velocity and temperature obtained from solving boundary layer equations for subsonic 
jet Mj=0.85 and Re=10 5 . The mean flow parameters are slowly varying with axial 
direction due to the high Reynolds number used in the simulation. The most unstable 
frequency and its subharmonics are obtained from the linear stability analysis for this 
mean flow as shown in figure (4). The disturbance growth rate is decreasing with axial 
direction due to the mean flow variations. The most unstable frequency calculated for 
this flow is around 1.5 which gives Strauhal number of 0.477 (st=f dj /Uj). 

A comparison between different McCormack type numerical schemes is done to decide 
which scheme is suitable for the noise prediction. Figure (5.a) represent a comparison 
between the linear stability theory , 4/2 and DRP/optimized schemes for too many points 
per wave length. The two schemes give good results for 17 number of points per wave 
length. As the amplitude of the disturbance becomes large, the results do not agree with 
the linear stability because the linear stability theory is no longer valid. By decreasing 
the number of points per wave length, some of these schemes is not performing well as 
shown in figure (5.b). For 9 points per wave length, the dispersion and dissipation errors 
are significant for the 4/2 scheme. Other schemes like 6/2, 6/4 and DRR/opt. still have 
good accuracy. 

Comparison between different outflow boundary conditions is introduced in figures (6) 
and (7). The reference solution is the long domain results where the computational 
domain is chosen to be long enough (x/rj=100) such that the disturbance becomes steady 
and periodic in the required domain (x/rj=60) before the disturbance hits the outflow 
boundary. It is clear from figure (7) that the perfectly matching layer technique gives 
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minimum reflections at the boundary and the buffer domain gives some reflections but it 
has less reflections than the characteristic boundary conditions proposed by Rudy 
(1980). 

The pressure disturbance contours are presented in fig. (8) for supersonic hot jet 
(Mj=2.5) with zero free stream velocity and temperature ratio equal 2.25. It is clear that 
the pressure disturbance propagates within a cone which means that the noise radiation 
from the instability waves is confined within a wedge while the noise radiation from 
subsonic jet has no specific pattern as clear from fig. (9). There is upstream influence of 
the pressure waves for subsonic jet and the jet radiates noise in every direction. The 
density and vorticity contours are shown in figures (10) and (11) where the shear layer 
roll-up and the vortex pairings are clear. 

The sound radiation by vortex roll-up and pairings for low Reynolds number subsonic 
jet (Re=2500, Mj=0.8) is investigated in figures (12)-(14). 

The jet is excited at the inflow by the most unstable frequency f and its subharmonics f/2 
and f/4. The Fourier transform of the energy amplitude is plotted in figure (12). It is 
clear from this graph that the second subharmonic is dominant . The energy grows in the 
linear region and then saturate at about x/rj=40 for the second subharmonic where the 
source of the noise is located. The amplitude of the disturbance starts to decay after the 
saturation where the jet becomes turbulent. The vorticity contours plotted in fig. (13) 
shows the vortex pairs where the amplitude of the second subharmonic becomes large. 
This means that the large vortical structures are responsible for the noise radiation. 

Using the wave equation, the sound pressure levels in the far field are calculated as 
shown in figure (14). No directivity pattern is predicted for the subsonic jet and the 
maximum sound pressure level obtained for a disturbance of amplitude 0.001 at the 
inflow is 140.5 decibels. 
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